function [rmse,stt,smooth_st,etamat,yfor,ynewmat]= gensys_sspace( Y , modstr , fullp , into , instateloc , outpath , namestr)
% =========================================================================
% =========================================================================
% GENSYS_SSPACE.m  
% 
% Plot the actual series vs. 1step ahead forecasts 
%      Real-time and smoothed states for the whole system 
%      innovations and table with analysis of their properties   
% 
% function gensys_sspace( Y , modstr , fullp , into , outpath , namestr)
%
% Inputs 
% -------
% FULLP         Vector of paramaters ( including calibrated ) 
% INTO          Structure with fields 
% INSTATELOC    Vector of Initial State Location 
% OUTPATH   Path to store graph and tables ( do not call or leave empty ) 
% NAMESTR   String appended to all files created here 
% 
% INTO.SERNAMES     Series names for Columns of Y 
% INTO.STNAMES   Names for the states, size NY  of G1 
% INTO.SHNAMES      Shock names for SDX 
% INTO.STARTEST     Start Date for Estimation 
% INTO.FINISHEST    Finishdate for Estimation 
% INTO.SAMPLE       Sample Vector 
% 
% Also use INTO to pass on all arguments to PLOTSAJ.M 
% 
% If OUTPATH and NAMESTR  provided, following files will be created 
% 
% INSTATELOC :  If empty then the initial state is set at zero 
%               Else, should indicate to which states should the first
%               observation in Y be matched 
%               E.g.  instateloc = [ 1 4 2 ] means the Y(1) will be matched
%               to these states and the filtering will start one
%               observation later 
%               Instatloc should be [1 nn] or [nn 1] 
%
% 
%
% Output 
% ======
% [YFORC namestr]      Plots of actual series vs. Forecasts 
% [INNOV namestr]      Plots of model innovations (SDX) 
% [ST namestr]         Plots of Real-time and smoothed states 
% [INNOV ANALYSIS namestr] Table with std, cross-correlations and VAR(1) 
%                           estimates for the innovations 
% 
% YNEWMAT                New Matrix which plots each series as the initial
%                        states are fixed to their smoothed values at time
%                        0 and the history of shocks are set to zero one at a time 
%                        allows to gauge the contribution of the shocks to
%                        the evolution of the endogenors variables. 
%                        [ T NX+1 NSER ] Last column in each page plots the values of T
%                        with all shocks "on" to see if there is an error
%                        in the computation. This usually leads to errors
%                        of 1e-10. 
%                       
%
% Alejandro Justiniano 9/7/2005 
% Added tables for the states and the innovations  9/26/2005 
% Plugges constants in states & forecasting 11/2/2005 
% Counterfactual with the shocks see GENSYS_COUNTER.M but ONLY if provide
% the initial state 
% Changed managament of the initial state which now must be provided from
% outside ( location ) as an extra input 
% October 10 2006 
% Compute RMSE for 1sha 
%         New plots for forecasts and tables TAB 
% ====================================================================
cucd = cd; 
[nobs,nn] = size(Y); 
[G1,CC,R,eu,SDX,Z]=feval(modstr,fullp);
nx=size(SDX,1);

if isequal(eu,[1;1])==0;error('Indeterminacy');end; 
if nargin < 6 | isempty( outpath ) == 1 ;
    flag_save = 0 ;
else
    flag_save = 1; 
end
% ==========================================
% 1. Check the settings for the initial state 
% =========================================
if isempty(instateloc) == 1 
    disp('Initial state will not be used'); 
    flag_instate=0 ;
else 
    instatloc=instateloc(:);  
    if length(instateloc) ~= nn 
        error('Dimension Mistmatch: INSTATELOC and columns of Y'); 
    end 
    flag_instate = 1; 
    if max( instateloc ) > size(G1,1); 
        error('INSTATELOC entry exceeds dimension of state vector'); 
    end
end 
% ================================================
% 2.Check settings for state,series and shocks names 
% =================================================
[e,into]=ch_field( into, 'sernames' ); 
if e==0 
    into.sernames = fnumcell( 'Series ',[1:nn] ); 
end 
[e,into]=ch_field( into, 'shnames' );  
if e==0 
    into.shnames = fnumcell( 'Shock ',[1:nx] ); 
end 
[e,into]=ch_field( into, 'stnames' ); 
temp=size(G1,1); 
if e==0 
    into.stnames = fnumcell( 'State ',[1:temp] );
end; clear temp; 
% ================================================ 
% 3. Settings begin and end estimation 
% ================================================
[e,into]=ch_field(into, 'sample' , ([1:nobs])' ); 
[e,into]=ch_field(into, 'startest' , 1 ); 
[e,into]=ch_field(into,'finishest', nobs ); 
[e,into]=ch_field(into,'titfont', 8); 

tbeg=find(into.sample==into.startest); 
tend=find(into.sample==into.finishest); 
% =============================================
% 4. Remove constant and filter the series 
% =============================================
flag_cons=0; 
if any( CC~=0 ) == 1 ;
    flag_cons=1; 
    Cz=Z*CC;
    Cz=repmat( Cz' , [size(Y,1) 1] ); 
    Y=Y-Cz;
end
[stt,etamat,smooth_st,yfor,vstar,logLnc]=...
    kfilter_full(Y,Z,G1,R,SDX,flag_instate,instateloc);

% =====================================================================
% Analyze contribution of each shock to the variability of each series 
% =====================================================================

% ===================================================================
% Truncate to the estimation sample or account for the initial state 
% ===================================================================
if tbeg==1 & flag_instate == 1 
    tbeg=2 ; 
end

% Add back the constants to the series and the states 
% ====================================================
if flag_cons == 1
    Y = Y + Cz;
    newt = size(stt,1);
    Cmat = repmat( CC' , [newt 1]);
    stt = stt + Cmat ;
    smooth_st = smooth_st + Cmat;
    yfor = yfor + Cz(1:newt,:);
    clear newt; 
end

% ===================================================================
%  Counterfactuals by shutting down innovations 
% ===================================================================
szero=smooth_st(tbeg,:); 
ezero=etamat(tbeg,:); 
if flag_instate== 1
    disp('Compute the counterfactuals'); 
    [ynewmat,difmat] = gensys_counter( Y(tbeg+1:end,:) , Z, G1, R, CC, etamat(tbeg+1:end,:) , szero );
    clear szero ezero
    if max( max( abs(difmat) ) ) > 1e-2;
        warning('GENSYS_COUNTER: Numerical problems do not match the generated Ys to the actual series');
    end
    tinto.legend={'Count','Actual'};
    tinto.cnames=tinto.legend;
    tinto.m = 3;
    tinto.n = ceil( nx / 3);
    tinto.pnames = into.shnames;
    tinto.xaxis = into.sample( tbeg+1:end );
    if flag_save == 1
        tinto.opath = outpath;
    end
    temp_t = size( ynewmat ,1 );
    ii=1;
    for ii=1:nn;
        tinto.title=[': ',into.sernames{ii}];
        count_mat = zeros(  temp_t , 2 , nx );
        count_mat(:,1,:) = reshape( ynewmat(:,1:nx,ii), [temp_t 1 nx] );
        count_mat(:,2,:) = repmat( Y( tbeg+1:end, ii) , [1 1 nx] );
        if flag_save == 1
            tinto.gname = ['Count ',into.sernames{ii}];
        end
        plotsaj( count_mat , 1 , tinto );
        disp(['Counterfactual for: ',into.sernames{ii}] );
    end
    clear tinto temp_t count_mat ii;
else 
    disp('WILL NOT DO the counterfactuals'); 
end 
% ================================
%  END counterfactuals 
% ===============================
if tbeg > 1; 
    stt=stt(tbeg:tend,:); 
    etamat=etamat(tbeg:tend,:); 
    yfor=yfor(tbeg:tend,:); 
    smooth_st=smooth_st(tbeg:tend,:); 
end 
samp=into.sample( tbeg:tend );
nobs=size( yfor ,1 ); 
Y=Y(tbeg:tend,:); 
if size(Y,1) ~= nobs 
    error('Dimension Mistmach Y and Forecast') 
end 

into.title = ' ,Actual & 1-SAF'; 
into.legend = {'Actual ';'1-SAF '}; 
into.cnames = into.legend;  
into.pnames = into.sernames; 
into.xaxis  = into.sample(tbeg:tend) ;


nsubs=3; 
hcont=zeros(nobs,nx,nn); 
temp=nsubs*ceil(nobs/nsubs); 
ipos=1:ceil(nobs/nsubs):temp;
ipos=[ipos nobs+1]; 

clear temp; 
cmat=[0 0 0;
    1 1 1;
    1 0 0;
    0 1 0.5;
    1 1 0.25;
    1 0.5 0;
    0 0 1;
    1 0.25 1
    127/255 1 212/255];

ZH=Z*R;
ii=1;
for ii=1:nn;
    temp=repmat(ZH(ii,:),[nobs 1]);
    hcont(:,:,ii)=temp.*etamat;
    figure;
    jj=1;
    for jj=1:nsubs;
        subplot((nsubs+1)/2,(nsubs+1)/2,jj);
        x1=into.xaxis(ipos(jj):ipos(jj+1)-1);
        y1=hcont(ipos(jj):ipos(jj+1)-1,:,ii);
        ypos=y1.*(y1>0); 
        yneg=y1.*(y1<=0); 
        bar(x1,ypos,'stacked'); %colormap(cmat(1:nx,:));
        hold on; 
        bar(x1,yneg,'stacked'); 
        %colormap(hot);
        xlim([x1(1)-0.25 x1(end)+0.25]);
        %dif=Y(ipos(jj)+1:ipos(jj+1)-1,ii)-Y(ipos(jj):ipos(jj+1)-2,ii);
        dif=Y(ipos(jj)+1:ipos(jj+1)-1,ii)-yfor(ipos(jj)+1:ipos(jj+1)-1,ii);
        plot(x1(2:end),dif);
        hold off;
    end;
    subplot((nsubs+1)/2,(nsubs+1)/2,nsubs+1);
    bar([1:2],zeros(2,nx),'stacked');
    legend(into.shnames,'Best');
    suptitle_withpatch([char(into.sernames(ii)),' Historical Contribution']);
    cd(outpath);
    print(gcf,'-deps',[char(into.sernames(ii)),' hist']);
    cd(cucd);
end;
% =================================================
%        Create output cells 
% ==================================================
samp_cell = num2cprec( samp , 2 ); 
temp_top = combine_cells( [] , ' ' , into.stnames , ' Stt ' )'; 
ocell   = crtcell( stt , samp_cell , temp_top , 'Stt States' , 1 ); 
temp_top = combine_cells( [] , ' ' , into.stnames , ' St Smooth ' )'; 
st_cell = crtcell( smooth_st , samp_cell  , temp_top , 'Smooth State' , 1 ); 
ocell   = merge_cells( ocell, st_cell , 2 ); 
if flag_save == 1 
savecell( ocell , cucd, outpath , ['STATES ',namestr] ); 
end 
clear ocell; 
temp_top = combine_cells( [] , ' ' , into.shnames , ' innov' ); 
ocell = crtcell( etamat , samp_cell , temp_top , 'Disturbance Innovation'); 
if flag_save == 1 
    savecell( ocell , cucd , outpath , 'Innovations' ); 
end 
% =================================================
% Plot the actual series relative to the forecasts 
% =================================================
if flag_save == 1 
    into.opath = outpath; 
    into.gname = ['YFORC ',namestr]; 
end 

if nn > 4     
    if nn < 6;
        into.m=2;
        into.n=3;
        into.titfont=8;
    else
        into.m=3;
        into.n=3;
        into.titfont=7;
    end 
end 

% =================================================
% Plot the actual series relative to the forecasts 
% =================================================
xmat = zeros( nobs, 2, nn); 
xmat(:,1,:) = Y ; 
xmat(:,2,:) = yfor ; 
into.width=[1.5 2]; 
into.style={'b-','r.-'}; 


rmse=xmat(:,1,:) - xmat(:,2,:); 
rmse=rmse.*rmse; 
rmse=sum(rmse,1); 
rmse=num2cprec(rmse); 
rmse=combine_cells([],' : ',into.sernames,rmse); 

set(0,'DefaultAxesFontSize',7);
count=0;
isub=0;
ii=1;
[junk,temptit]=extrsubf(outpath,cucd); 

forcell=emptycell(length(into.xaxis)+1,2*nn+1); 
forcell(1,1:2)={'Actual','Forecasts'}; 
forcell(2:end,1)=num2cprec(into.xaxis); 

for ii=1:nn+1;
    if isub==0;
        figure;
    end;
    isub=isub+1; 
    if ii <= nn
        subplot(into.m,into.n,ii);
        plot(into.xaxis,squeeze(xmat(:,1,ii)),'c-','LineWidth',1.5);
        hold;
        plot(into.xaxis,squeeze(xmat(:,2,ii)),'k-','LineWidth',2);
        hold off;
        title(rmse{ii},'FontSize',into.titfont);
        xlim([into.xaxis(1) into.xaxis(end)]); 
        forcell(2:end,ii*2:ii*2+1)=num2cprec(xmat(:,1:2,ii)); 
        forcell(1,ii*2)={[into.sernames{ii},' Actual']};  
    else
        subplot(into.m,into.n,ii);
        plot(zeros(10,nn));
        lh=legend(rmse,'Location','Best' ) ;
        set(lh,'FontSize',7);
    end
    if isub==(into.m*into.n) | ii==(nn+1);
        isub=1;
        count=count+1;
        cd(outpath);
        if into.m*into.n <= nn+1
            tempname='YFORC';
        else
            tempname=['YFORC ',num2str(count)];
        end
        suptitle_withpatch(['1SAF (dark) and RMSE',temptit]); 
        saveas(gcf,tempname,'fig');
        print(gcf,'-depsc',tempname);
        cd(cucd);
    end
end
set(0,'DefaultAxesFontSize','remove');

forcell=merge_cells(rmse,forcell,1); 
savecell(forcell,[],outpath,'Tab 1SAF'); 

clear ii isub lh count tempname junk temptit forcell tempcell;

clear xmat ; 
into=rmfield(into,'width'); 
into=rmfield(into,'style'); 

if flag_save == 1 
    into.gname = ['INNOV ',namestr]; 
end 


temp2=size(SDX,1); 
if temp2> 4     
    if temp2 < 6;
        into.m=2;
        into.n=3;
        into.titfont=8;
    else
        into.m=3;
        into.n=3;
        into.titfont=7;
    end 
end 
into.title=' ,INOV '; 
into=rmfield(into,'legend'); 
into=rmfield(into,'pnames'); 
into.cnames=into.shnames ; 
plotsaj(etamat,2,into); 

% =================================================
% Plot the actual series relative to the forecasts 
% =================================================
temp = size( G1 ,1 ); 
if temp > 8 
    into.m = 2; 
    into.n = 3; 
end
% =================================================
% Plot the actual vs. Smoothed States 
% =================================================
into.pnames=into.stnames; 
if flag_save==1 
    into.gname=['ST ',namestr]; 
end 
xmat=zeros(nobs,2,temp); 
xmat(:,1,:)=stt; 
xmat(:,2,:)=smooth_st;

into.legend = {'Stt ','Smooth'}; 
into.cnames = into.legend; 
into.title = '';
into.width=[1.3 1]; 
plotsaj( xmat , 1 , into ); 
clear xmat temp; 

% ===============================
%  Analysis of the innovations 
% ===============================
emat_cor=corrcoef(etamat); 
emat_std=std(etamat); 
[emat_var,emat_res]=ols(etamat(2:end,:),etamat(1:end-1,:));

emat_cell = emptycell( 4 , nx +1 ); 
emat_cell(1,1:2) ={'SAMPLED','INNOVATIONS'}; 

emat_cell(3,1) = {'STD INNOVATIONS'}; 
emat_cell(3,2:end) = ( into.shnames(:) )'; 
emat_cell(4,2:end) = num2cprec( emat_std ); 

emat_c2 = crtcell( emat_cor , into.shnames , into.shnames , 'CORR INNOV' , 1); 
emat_cell = merge_cells( emat_cell , emat_c2 , 1 ); 
clear emat_c2; 

emat_c3 = crtcell( emat_var , into.shnames , into.shnames , 'VAR1 INNOV' , 1 ); 
emat_cell = merge_cells( emat_cell , emat_c3 , 1 ); 
clear emat_c3; 

printcell( emat_cell ); 
if flag_save == 1 ; 
    savecell( emat_cell, cucd, outpath , ['Tab INNOV ',namestr] ); 
end